Read the data

df

Variables of interest

  • Tmax: max temperature
  • Tmin: min temperature
  • Tmed: mean temperature
  • Rain: percepitation
  • Sm : soil moisture
  • Iaf : leaf area index

- Tmed: mean temperature

bins = 9
Year = 1041
v = 'Tmed'
 
plot_density_(df = df, Year = Year,  v = v , bins = bins)
## Picking joint bandwidth of 1.05

plot_ts_meteo_(df = df,  v = v, scaleLimits = 5) 

- Rain: Total rain each year

bins = 9
Year = 1041
v = 'Rain'



 p <-
    df %>%
    select(var, eval(v), Year, var) %>%
    rename(region = var,
           value = eval(v)) %>%
    group_by(Year, region) %>% 
    #summarise(value = mean(value, na.rm = TRUE), .groups = 'drop') %>% 
    summarise(value = sum(value, na.rm = TRUE), .groups = 'drop') %>% 
    
    
    ggplot(aes(x = Year, y = value, fill = region, color = region)) +
    geom_line( ) +
    expand_limits(y = 0 ) +
    ylab(v) + xlab('year') +
    theme_ipsum() +
    labs(title = paste('Total ',v,' per year in both regions'))
  
  # Turn it interactive with ggplotly
  ggplotly(p)

- Sm:Soil moisture

bins = 9
Year = 1041
v = 'Sm'
 
plot_density_(df = df, Year = Year,  v = v , bins = bins)
## Picking joint bandwidth of 20.6

plot_ts_meteo_(df = df,  v = v, scaleLimits = 5) 
v = 'Sm'

data_xts <- 
    df %>%
    select(var, eval(v), Year, var,Date ) %>%
    mutate(time_index = as.yearmon(Date)) %>% # create a time index for the data
  
   # mutate(Month = month.name[Month]) %>% 
    rename(region = var,
           value = eval(v)) %>%
    group_by(Year, region, time_index) %>% 
    #summarise(value = mean(value, na.rm = TRUE), .groups = 'drop') %>% 
    summarise(value = median(value, na.rm = TRUE), .groups = 'drop') %>% 
    
    spread(key = region, value = value)

ggseasonplot(
  ts(data_xts$rdd, start = c(min(data_xts$Year)), frequency = 12),
  polar = TRUE,
  year.labels = FALSE,
  year.labels.left = FALSE
) +
  ylab("Soil moisture") +
  ggtitle("Polar seasonal plot: Soil moisture - RDD region") +
  theme(legend.position = "none")

ggseasonplot(
  ts(data_xts$rvv, start = c(min(data_xts$Year)), frequency = 12),
  polar = TRUE,
  year.labels = FALSE,
  year.labels.left = FALSE
) +
  ylab("Soil moisture") +
  ggtitle("Polar seasonal plot: Soil moisture - RVV region") +
  theme(legend.position = "none")
## Warning: Removed 104 rows containing missing values (`geom_line()`).

  df_sm <- 
  df %>%
    select(var, eval(v), Year, var, Month) %>%
    rename(region = var,
           value = eval(v))

 df_sm %>% 
    #mutate(Month = as.yearmon(Date)) %>% # create a time index for the data
    group_by(Year, region) %>% 
    summarise(value = mean(value, na.rm = TRUE), .groups = 'drop') %>% 
    ggplot(aes(x = Year, y = value, fill = region, color = region)) +
    geom_line( ) +
    expand_limits(y = 0 ) +
    ylab(v) + xlab('year') +
    theme_ipsum() +
    labs(title = paste('Mean ',v,' per year in both regions')) +
  facet_wrap(~region) 

 # df_sm %>% 
 #   filter(region == 'rdd') %>%
 #   group_by(Month, region, Year) %>%
 #   summarise(value = mean(value), .groups = 'drop') %>%
 #   ggseasonplot(ts(value, start = c(min(Year)),frequency = 12),
 #                year.labels = TRUE,
 #                year.labels.left = TRUE) +
 #   ylab("$ million") +
 #   ggtitle("Seasonal plot: antidiabetic drug sales")

  # Turn it interactive with ggplotly
  #ggplotly(p)
  • Iaf : leaf area index

Prepare the data

# ## The xts() function creates a time series object. The order.by parameter uses the dates from the data, which are converted to R Date objects by the as.Date() function.
# 
# historical_rdd = xts(df_rdd[,c("Tmed", 'Rain')], order.by=as.Date(df_rdd$Year, df_rdd$Month ))
# 
# plot(historical_rdd$Tmed)
# 
# ## The ts_regular() function gives the time series a regular (daily) interval by adding NA values for missing dates. A regular interval will be needed later for decomposition analysis.
# 
# historical_rdd = ts_regular(historical_rdd)
# 
# 
# ## The na.fill() function fills in those missing dates by duplicitating (extending) values from previous days.
# 
# historical_rdd = na.fill(historical_rdd, "extend")
# 
# ## The window() function clips off the starting and ending dates so the number of years covered is a multiple of four. This will be needed later when the data needs to be aggregated into monthly periods.
# 
# historical_rdd = window(historical_rdd, start=as.Date("2000-01-01"), end=as.Date("2020-12-31"))